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The stationary points (SPs) of the potential energy landscapes (PELs) of multivariate random po¬ 
tentials (RPs) have found many applications in many areas of Physics, Chemistry and Mathematical 
Biology. However, there are few reliable methods available which can find all the SPs accurately. 
Hence, one has to rely on indirect methods such as Random Matrix theory. With a combination 
of the numerical polynomial homotopy continuation method and a certification method, we obtain 
all the certified SPs of the most general polynomial RP for each sample chosen from the Gaussian 
distribution with mean 0 and variance 1. While obtaining many novel results for the finite size case 
of the RP, we also discuss the implications of our results on mathematics of random systems and 
string theory landscapes. 


I. INTRODUCTION 

The surface drawn by a potential, V(x), with x = 
(xi,..., xat), is called the potential energy landscape (PEL) 
of the corresponding physical or chemical system [1] [5]. The 
stationary points (SPs) of a PEL, defined by the solutions of 
the equations dV{x)/dxi = 0 for i = 1, ... ,N , provide im¬ 
portant information about the physics and chemistry of the 
corresponding system. If the parameters or the coefficients 
of V(x) are chosen from some random distribution, then 
V(x) is called a random potential (RP). In recent years 
more attempts to understand the statistical patterns be¬ 
hind SPs for random PELs have been made because of its 
applications in such diverse areas as string theory II ll]> 
cosmology (see e.g. , |5HI3), statistical mechanics HMZ], 
neural networks m- Similar problems appear in statistics 
|19| and in a pure mathematics context, e.g. in topology 
|20H22] . In fact, the mathematical question how many real 
solutions, on an average, does a random system of poly¬ 
nomial equations have is a classic problem. For random 
univariate polynomials of degree D with coefficients taking 
i. i. d. values from the Gaussian distribution with mean 0 
and variance 1, called the Kac formulation, the mean num¬ 
ber of real solutions is f In U for D ^ oo |23l . If the 

f-th coefficient of a random polynomial is allowed to take 
i. i. d. values from the Gaussian distribution with mean 0 
and variance (^), called the Kostlan formulation, then the 
mean number of real solutions is VD as D ^ oo [25] . 
There have been various attempts to extend these results 
to the multivariate case |27H3l]. 

However, progress in studying the statistics of the SPs of 
PELs has been slow due to conceptual as well as techni¬ 
cal problems. First, from a computational viewpoint, the 
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stationary equations are usually nonlinear. Hence, solving 
these equations is usually extremely difficult. One may use 
the Newton-Raphson method, its sophisticated variants, or 
many other numerical methods to solve these equations. 
These methods can find many SPs at best, as opposed to 
all of them. Theoretically, the most recent progress was 
based on successful applications of random matrix theory 
(RMT) for exploring the PELs of certain types of RPs. In 
this approach, the Hessian matrix Hij = d^V(x)jdxidxj, 
is treated as a random matrix. Then, standard RMT results 
can be employed to extract valuable information about 
the eigenvalue distribution of H. In recent years, there 
are many results available on the probability that random 
matrices of various types has indefinite spectrum, see e.g. 
|32I[M] . However, constraining this analysis for the SPs 
rather than on arbitrary points of the iV-dimensional space 
amounts to taking the stationary equations under consid¬ 
eration into the matrix approach. This can be successfully 
done only under assumptions of gaussianity, and either con¬ 
ditions of statistical isotropy and translational invariance 
[111 [121 mi (though this condition can be made milder in 
certain situations, see e.g. , Sec. 4 in [3S]) or isotropy and 
restriction to a spherical surface [T51 [T51 [55] . The latter 
context is natural for studying various aspects of glassy 
transitions such as the number of minima and their com¬ 
plexity, mainly in the thermodynamic limit of high dimen¬ 
sion as TV —)■ 00 where RMT produces the most explicit and 
universal results. 

There are, however, a few issues with this approach as well. 
Many cases exist where N is actually a fundamental param¬ 
eter in the physical description and may be finite. Indeed, 
many physical, chemical or biological systems have a finite 
number of fields, particles, neurons, etc, which is repre¬ 
sented as the number of variables N. Moreover, with the 
RMT approach, it is not yet possible to obtain further in¬ 
formation about an individual SP. Gomputing the variance 
of the number of SPs using RMT is also quite difficult, and, 
as of yet, is still largely an unsolved problem |55|. We note 
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that recent progress has been made in which the average 
number of SPs of each index (defined as the number of 
negative eigenvalues of the Hessian matrix evaluated at the 
given SP), for any finite N, is computed analytically |14f - 
[IZ1I3S]. The shortcoming of these approaches is that the 
theory only works for the gaussian isotropic RPs with either 
spherical constraints or statistical translational invariance. 
In this article, we provide a numerical scheme to overcome 
all the difficulties to find the SPs of a completely generic 
RP. We use the numerical polynomial homotopy contin¬ 
uation method which finds all the isolated solutions of a 
multivariate system of polynomial equations with proba¬ 
bility one |37I[39| . To strengthen our results, we also use a 
certification approach which certifies if a given numerical 
approximate corresponds to an exact distinct root of the 
system. The certification approach is known as Smale’s a- 
theory and certifies that the numerical approximation is in 
the quadratic convergence region of a unique nonsingular 
solution of the system |40fH4] . 

In the remainder of the paper, after specifying our RP, we 
describe the computational methods used in our work in 
Section]^ In Section [HT] we present our results and discuss 
their implications. Finally, in Section |IV| we conclude. 

II. THE MOST GENERAL POLYNOMIAL 
RANDOM POTENTIAL AND NUMERICAL SET 

UP 

The most general random multivariate potential with poly¬ 
nomial nonlinearity is 

V{x)= ^ aaX'^\..x%^ (1) 

\a\<D 

where N is the number of scalar fields, D is the highest 
degree of the monomials, a = (oi,..., a^) G is a multi¬ 
integer, and |a| = ai + ... +a]v- are random coefficients 
which take i. i. d. values from the Gaussian distribution 
with mean 0 and variance 1, known as the Kac formulation 
for the multivariate case. 

Finding all the SPs of V{x) boils down to simultaneously 
solving the system of equations = 0, i = 1,...,1V, 

each of which is a polynomial equation of degree D—\. The 
classical Bezout theorem states that for generically chosen 
coefficients, the number of complex (which include real) 
SPs, counting multiplicities, is the product of the degrees 
of each of the stationary equations, {D — 1)^. Using this 
fact, we employ the numerical polynomial homotopy con¬ 
tinuation (NPHC) method which guarantees that we will 
find all the complex numerical approximates of a system of 
multivariate polynomial (or, having polynomial-like non¬ 
linearity) equations |37f - l5^ . The NPHC method has been 
used to explore potential energy landscapes of various prob¬ 
lems arising in physics and chemistry gSHSl]. Here, one 
first creates an easier to solve system which has the same 
variables as the original system as well as the same num¬ 
ber of solutions as the classical Bezout count; then, each 
solution of the new system is tracked to the original system 
with a single parameter to obtain all the numerical approxi¬ 
mates of the original system. It is rigorously proven |55[I56| 


that by tracking the paths over complex space one can be 
guaranteed to find all the isolated complex solutions of the 
system using the NPHC method. The interested reader is 
referred to the early references I37HM1I15H511 for a detailed 
description of the method. For the NPHC method, we use 
the Bertini software 123 EH- 

Though the NPHC method guarantees one will find all nu¬ 
merical solutions, skeptics may wonder if the numerical so¬ 
lutions are good enough. Smale and others defined a “good 
enough” numerical solution of a system of polynomial equa¬ 
tions as a point which is in the quadratic convergence region 
of a nearby exact solution of the system 123 mi- Once the 
numerical solution is good enough in the above sense, it can 
then be approximated to arbitrary accuracy. Such a numer¬ 
ical solution is called a certified solution. Smale also showed 
how to certify a numerical solution using only data such as 
the Jacobian and higher derivatives at the numerical solu¬ 
tion via a-theory. A nontrivial result of Smale is that if a, 
which is computed using the Jacobian and higher deriva¬ 
tives of the system of equations, at a given numerical solu¬ 
tion is less than (13 — 3-\/7)/4, then the numerical solution 
is within the quadratic convergence region of the nearby 
exact solution of the system. Hence, this approximate is a 
certified numerical solution. We use an implementation of 
a-theory provided in alphaCertified |42| in the present 
work (see Refs. 03122] for certification of the PELs). We 
call a solution a real finite solution if the imaginary part of 
each variable at the solution is 10“^°, the maximum eigen¬ 
value of the Hessian matrix at that solution is < 10^^, and 
it is certified using alphaCertif ied. 

In practice, we must restrict ourselves to N = 2,..., 11 and 
D = 3,... ,5 and the sample size for each pair of {N, D) is 
1000. The limitation comes due to the combined computa¬ 
tion of solving equations, certifying solutions and comput¬ 
ing Hessian eigenvalues. 


III. RESULTS AND DISCUSSION 

In this Section, we present our results from our numerical 
experiments. We first find all the isolated complex SPs 
of V{x) for each of the 1000 samples, and then compute 
various quantities from the SPs. We present the results for 
Y = 2 ... 11 and D = 3.. .5. We also discuss the results in 
a physical and mathematical context. 


A. Average number of real SPs 

The number of complex solutions of the stationary equa¬ 
tions for generic coefficients is always {D — 1)^ for our RP. 
A more interesting quantity is the mean number of real 
SPs which we plot as a function of N for different values 
of D in Figured For limited values of U, a similar plot 
was drawn in |58] for a different RP which was a Taylor ex¬ 
pansion and coefficients drawn from a uniform distribution. 
Remarkably, we see similar behaviour of the mean number 
of SPs, namely, the number grows rapidly as N increases. 
In fact, in Figure [T] we compare our numerical results for 
the mean number of real SPs as a function of N with an 
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analytically computed upper bound, V^{D — 
computed in m where the coefficients were i. i. d. picked 
from the Gaussian distribution with mean 0 and variance 
being a multinomial, called Kostlan-Shub-Smale formula¬ 
tion I1SJI27]. While the numerical results are in fact far 
below from the upper bound, this should be expected since 
it is known that mean number of real SPs for the Kostlan 
formulation is significantly higher than the Kac formulation 

m- 

As is true with most other analytical computations related 
to random systems, analytically computing the variance of 
the number of SPs is almost a prohibitively difficult task. 
Figure shows that the variance also behaves in qualita¬ 
tively the same fashion as the mean indicating that as N in¬ 
creases the mean of the number of real SPs are very spread 
out from each other and from the mean. To further inves¬ 
tigate the spread, we plot histograms of the number of real 
SPs in Figures [3]|^ which show that the peak in the his¬ 
tograms is indeed around the mean value of the number of 
SPs, i.e., a unimodal distribution. Away from the peak, the 
histograms spread out smoothly. However, our results sug¬ 
gest that these histograms are not always symmetric with 
respect to the peak but often exhibit right-skewedness. 
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Figure 3. Histogram of the numbers of SPs for D = 3. A = 5, 6 
in the top row, A = 7, 8 in the middle row, and A = 9,10 in 
the bottom row. 



Figure 1. Mean number of SPs as a function of A for various val¬ 
ues of D. Here, ’UB’means the upper bound -\/2(F) — 

m- 


Variance of No. real SPs 
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Figure 2. Variance of the number of SPs as a function of A for 
various values of D 




Figure 4. Histogram of the numbers of SPs for D = 4. A = 5, 6 
in the top row, A = 7, 8 in the middle row, and A = 9,10 in 
the bottom row. 


B. Average number of minima 

Since the potential may be unbounded from below (or 
above), we must clarify that by minima we simply mean 
SPs at which all the Hessian eigenvalues are positive 
definite. Though the mean number of real SPs in¬ 
creases exponentially with increasing A, the mean num¬ 
ber of minima decays exponentially with increasing N, 


as shown in Figures and We compare our results 
with an analytically computed upper bound |59| . namely, 

( _ /vr2 In 3 i (-^ + 1) l v, f r)_i 

ATel 4+2 f ^ where AT is a positive constant. 
Since a prescription for computing K is not available in 
[59], we obtain it by fitting the upper bound with the nu¬ 
merical data. Our results obeys the upper bound. In fact, 
the mean number of minima itself appears to qualitatively 

behave as = in(D-i)) ^ must point out that 

as a function of D, while keeping N fixed, the mean num- 
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Histogram of «(SPs). N=3 D=5 


Histogram of «(SPs), N=4 D=5 
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Figure 5. Histogram of the numbers of SPs for D = 5. N = 3,4 
in the top row, and = 5, 6 in the bottom row. 


ber of minima appears to be a slowly increasing function. 
This is not surprising since in this case the total number 
of complex solutions themselves {D — 1)^ increases dras¬ 
tically. We do not have enough data in D to extract the 
precise behaviour of the increase. However, the aforemen¬ 
tioned upper bound with now fixing N and varying D im¬ 
plies that the increase of mean number of minima should 
be much slower than the increase when D is fixed and N is 
varied because of the natural logarithm term. The variance 
of the number of minima also exponentially decays as N in¬ 
creases indicating that the mean number of minima gives 
a good estimate for the actual number of minima for any 
randomly picked sample. We anticipate that our numerical 
results will instigate an anlytical study for the variance of 
the number of real SPs and minima. 


Mean No. of Minima 
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Figure 6. Mean number of minima and comparing an analyti¬ 
cally computed upper bound m on the mean number of min- 


Variance of No of Minima 



D = 3 
D = 4 
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Figure 7. Variance of the number of minima. Since the variance 
is 0 for N > 7, we plot the data in the regular scale rather than 
the log scale. 


C. Average number of minima at which V > 0 (i.e., 
the de Sitter minima) 


The number of minima at which V > 0 may be viewed as 
de Sitter vacua if the potential is considered to mimic the 
statistical aspects of the string theory landscapes. Figures 
and 1^ essentially exhibit the same qualitative behaviour 
as the number of minima. The implications of these results 
on the string theory landscapes may be important, i.e. , the 
number of de Sitter vacua decreases exponentially with in¬ 
creasing the dimension of the moduli space, which is similar 
to the conclusion laid out in |5S] (see also |71[S]): the num¬ 
ber of vacua of a RP with tunneling rates low enough to 
maintain metastability exponentially decreases as a func¬ 
tion of the moduli space dimension. Though in the latter 
work, the RP was a bit different than ours in that in the 
former the RP was a Taylor expansion of a generic func¬ 
tion, with random coefficients i. i. d. picked from uniform 
distributions from specific ranges, remarkably, the overall 
conclusion remains the same even with our investigations 
with more general type of vacua (de Sitter vacua without 
any further conditions on bubble nucleation rate or else). 
However, we also point out that as we increase the degree 
of the potential (which would correspond to the truncation 
degree in the Taylor expanded version of |58| ) the number 
of de Sitter vacua does rise exponentially though with a 
linear rate. 

Mean No of de Sitter minima 
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Figure 8. Mean number of de Sitter vacua 


Variance of No of de Sitter minima 
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Figure 9. Variances of the number of de Sitter vacua 


D. Average number of Index-sorted SPs 

In Figures and El we further explore the properties of 
the real SPs by sorting them according to their Hessian in¬ 
dices. Since for the RP, generic random samples do not pos¬ 
sess singular solutions according to Sard’s theorem [60], the 
number of negative eigenvalues of the Hessian is indeed the 
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correct Morse index of the given real solution. Because we 
have moderate size systems in N and D at our disposal, we 
can not predict the precise behaviour for general N and D. 
However, it is clear from our results that the plots for I vs 
number of real SPs of index I tend to be bell-shaped curves. 
Such a plot is observed for many other potentials such as 
the Lennard-Jones potential |61l[62] . nearest-neighbour XY 
model spherical 3-spin model jM], the Ku- 

ramoto model |SS], etc. In jHl], such a behaviour of the 
number of real SPs with a given index for general potential 
was analytically derived. In short, in the landscape of RPs, 
there are vast number of SPs with non-zero index compared 
to the number of minima, which has also been analytically 
observed in random set ups |59| . 



Figure 10. Mean number of SPs as a function of I/N, where 
I is the Hessian index which runs from 0 to N, for D = 5 and 
various N. The corresponding plots for other values of D exhibit 
the same qualitative behaviour. 


not likely to have zero eigenvalues, though there are indeed 
SPs having eigenvalues closer to zero, creating the cleft at 
the zero eigenvalue in the histograms. From the statistical 
physics point of view, as argued in [7], a possible reason 
of the cleft may be that the Hessian matrices may have a 
component that can be described by the Altland-Zirnbauer 
Cl ensemble |66], though investigating this aspect further 
is beyond the scope of the present work. 
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Histogram of all Hessian eigenvalues (-40,+40). N=7 D=3 


8000 

6000 

4000 

2000 


12000 

10000 

8000 

6000 

4000 

2000 


-20 0 20 40 

Histogram of all Hessian eigenvalues (-40,+40), N=8 D=3 
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Histogram of all Hessian eigenvalues (-40,+40), N=9 D=3 
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Figure 11. Variance of the number of SPs of index 7 as a function 
of I/N, where I is the Hessian index which runs from 0 to N, 
for D = 5 and various N. The corresponding plots for other 
values of D exhibit the same qualitative behaviour. 


E. Histograms of the Hessian eigenvalues 


In Figures p2|T4 we plot histograms of the Hessian eigen¬ 
values where the x-ranges in the plots are chopped beyond 
±40 to show the behaviour of the plots in the middle range. 
Such histograms of eigenvalues of RPs are widely studied, 
usually by considering the Hessian matrices of RPs as ran¬ 
dom symmetric matrices and then applying Wigner’s semi¬ 
circle law or other results. Our results, however, are distinct 
in that they represent the histograms of Hessian eigenvalues 
computed at all the SPs for each sample. The clefts in the 
histograms are a departure from Wigner’s semicircle law. 
The reason behind the cleft can be explained using Sard’s 
theorem m, which yields that our dense random system 
of stationary equations will not possess singular solutions 
for almost all values of coefficients. This means that we are 


Figure 12. Histogram of eigenvalues, D — 3 
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Hblogram of iiJl Hessian eigenvalues (-40,+40), N=5 D=4 
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Figure 13. Histogram of eigenvalues, U = 4 


F. Histograms of the lowest Hessian eigenvalues 


The lowest eigenvalue of the Hessian matrix evaluated at a 
real SP is one of the most important physical quantities as it 
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Figure 14. Histogram of eigenvalues, D = 5 


determines physics up to certain extent. Many interesting 
properties of the lowest eigenvalues of different potentials 
also have deep connection to catastrophe theory |57]. In 
Figures [TsflTl we plot histograms of the lowest eigenvalues 
computed at all real SPs for all samples for various values 
of N and D. Our results yield that the tail of the lowest 
eigenvalues is long, though there are only a few events at 
the far end. 
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Histogram of the lowest eigenvalues D = 3. N = 9 
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Figure 15. Histogram of the lowest eigenvalues for D — 3. Data 
for smaller eigenvalues than the shown in the histograms are 
chopped away for the presentational purposes. 


G. Real vs Imaginary Plots of V at Complex 
Solutions 


In most applications we require the knowledge of real SPs, 
however, motivated by |53L I68| . we plot real vs imaginary 
parts of the potential at complex SPs in Figures for 
0 = 4, which yield that there are many complex SPs at 
which the potential V itself is a real quantity. Though 
in many physical models, both coordinates and potential 
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Figure 16. Histogram of the lowest eigenvalues for D = 4. Data 
for smaller eigenvalues than the shown in the histograms are 
chopped away for the presentational purposes. 
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Figure 17. Histogram of the lowest eigenvalues for D = 5. Data 
for smaller eigenvalues than the shown in the histograms are 
chopped away for the presentational purposes. 


have to be strictly real quantities and the complex solutions 
may appear just due to the complexification of the system 
in order to be able to use the complex algebraic geometry 
methods, in string phenomenology complex moduli do have 
physical meaning. We hope that our observations put for¬ 
ward in these figures may provide interesting avenues for 
further research to understand and interepret complexified 
fields. 


H. Average Timing 


Although not directly attacked, our experiments are re¬ 
lated to the 17th problem from the list of the eighteen un¬ 
solved problems for the 21st century that Smale laid down 
in Ref. [69] which we quote verbatim: “Can a zero of N 
complex polynomial equations in N unknowns be found 
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Im(V) against Re(V).N=5 D=4 


Im(V) against Re(V),N=6 D=4 
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Figure 18. Real vs Imaginary values of V at all complex solu¬ 
tions, D = 4 


approximately, on the average, in polynomial time with a 
uniform algorithm?” 

In order to gain some insight into this difficult problem, we 
consider the average time the Bertini 1.4 software takes to 
find one solution and consider this value for fixed fixed D. 


For Z? = 3 — 5, we note that for the systems that we con¬ 
sidered, the average time to track each path seems to grow 
exponentially in N as seen in Figure [T^ However, it is very 
important to note that the timings here includes the run¬ 
time of the homotopy path-tracking algorithm as well other 
nonessential functions. In addition, the timing information 
in Figure 19 was obtained using a single core processor (a 
2300 MhZ AMD processor) to avoid timing involved during 
parallel processes, though all the results presented in the 
previous subsections were carried out using a computing 
cluster of 64 processor. Hence, Figure [T9| in no way gives a 
good estimate of the actual averate time for path tracking, 
but just gives a very crude estimate in the hope to provide 
a rough guide for the practitioners. 


t (seconds ) 




Figure 19. N vs average time (in seconds) to track each solution 
path on a single machine. As explained in the text, this plot 
only provides a very crude estimate on average timing. Here, 
we start from A = 1 to get more data points. 


IV. CONCLUSIONS AND OUTLOOK 


In this article, we find all the critical points, complex and 
real, of the most general polynomial potential whose coef¬ 
ficients take i. i. d. values from the Gaussian distribution 
with mean 0 and variance 1. With this set up, we extract 
statistical results on random potentials (RPs). We have 
employed the numerical polynomial homotopy continuation 
(NPHC) method which guarantees that we have found all 
the critical points. Since the NPHC is a numerical method, 
choosing specific values for the tolerances to classify when 
the solutions are real becomes a delicate issue when TV and 
D are large enough. Hence, we combine a certification pro¬ 
cedure which provides a sufficient condition as to when a 
given approximate solution converges to a distinct real fi¬ 
nite solution using Smale’s a-theorem |4fil 1^ . Together 
with this fact and that the stationary equations of our RP 
always have (D — l)^ complex critical points for generically 
chosen coefficients, our numerical results can be treated as 
good as exact results. 

Due to various involved computations (solving equations 
using the NPHC method, certifying and sorting them, and 
computing Hessian eigenvalues), and limited computational 
resources, we have restricted ourselves to the number of 
variables N G {2,..., 11} and degree D G {3,..., 5}, and 
1000 samples for each (A, D). 

We have shown that the mean number of real SPs expo¬ 
nentially increases, whereas the mean number of minima 
exponentially decreases, with increasing N for various val¬ 
ues of D. Our results also yield that the mean number of 
minima increases for fixed N while increasing D, though 
the increase is linear. With an additional constraint on 
the minima that the potential evaluated at minima is pos¬ 
itive definite, we investigate the statistical aspect of the 
string theory landscape where our results may have im¬ 
portant consequences in counting the string vacua. These 
results compare well with an analytically computed upper 
bounds for these two average quantities j^. One of the 
main achievements of this paper is a numerical computa¬ 
tion of the variance of the number of real SPs and minima 
since it is prohibitively difficult to obtain analytical results 
for the variance. We observed that the variances exhibit 
the same qualitative behaviour as the mean. To further 
investigate the spread of the number of real SPs, we also 
plotted histograms of the number of real SPs which exhibit 
clear unimodality and right-skewedness. 

By sorting the real SPs according to their Hessian indices, 
we show that the mean number of real SPs having index 
I -I- 1 is significantly more than the those having index I. 
Extrapolating our results to higher N and D yields that 
the plot of I/N vs mean number of SPs with index I tend 
to be a bell-shaped curve, which appears to be a com¬ 
mon feature for many other potential energy landscapes 
such as the Lennard-Jones clusters HD ED, XY model 
gsi ig Ea [631 [65], spherical 3-spin mean-field model |64[ . 
etc. A milder version of this result, that the ratio between 
the number of SPs with index / > 0 and the number of min¬ 
ima exponentially blows up when N increases, has gained 
attention recently in both mathematics and string theory 

iziEn]. 


























Yet another novel result of our investigation is histogram of 
the Hessian eigenvalues computed at all real SPs since it is 
a departure from the traditional eigenvalue histogram stud¬ 
ies which consider Hessian matrices evaluated at arbitrary 
points and may arrive at Wigner’s semicircle law or else. 
Our histograms show clear bimodal symmetric behaviour 
with the cleft near the zero eigenvalue is pronounced as N 
increases. We explained this result using Sard’s theorem. 
We also observe that the histograms of lowest eigenvalues 
are long-tailed. 

Though usually physically interesting solutions are the real 
solutions, we observe in our numerical experiments that 
the potential may be real when evaluated at many of the 
complex (i.e., non-real) solutions. Such a phenomenon is 
observed and discussed previously for different potentials 
m 168] . Since in string theory, complex moduli play a 
prominent role in determining the physics of the model, 
our results may prompt a different interpretation of the 
feasible SPs. 

We also provided average timing information for our com¬ 
putation for finding solutions using homotopy continuation 
method in Figure |19| Though we do not claim to have di¬ 


rectly worked on Smale’s 17th problem |69| . nor to have an 
accurate timing estimate due to the overlaps with several 
other computations we anticipate that our crude estimates 
will guide future computations of the other practitioners of 
the homotopy continuation methods. 
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